home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The Very Best of Atari Inside
/
The Very Best of Atari Inside 1.iso
/
mint
/
mntlib25
/
atof.c
< prev
next >
Wrap
C/C++ Source or Header
|
1992-10-01
|
21KB
|
851 lines
/*
* double strtod (str, endptr);
* const char *str;
* char **endptr;
* if !NULL, on return, points to char in str where conv. stopped
*
* double atof (str)
* const char *str;
*
* recognizes:
[spaces] [sign] digits [ [.] [ [digits] [ [e|E|d|D] [space|sign] [int][F|L]]]]
*
* returns:
* the number
* on overflow: HUGE_VAL and errno = ERANGE
* on underflow: -HUGE_VAL and errno = ERANGE
*
* bugs:
* naive over/underflow detection
*
* ++jrb bammi@dsrgsun.ces.cwru.edu
*
* 07/29/89:
* ok, before you beat me over the head, here is a hopefully
* much improved one, with proper regard for rounding/precision etc
* (had to read up on everything i did'nt want to know in K. V2)
* the old naive coding is at the end bracketed by ifdef __OLD__.
* modeled after peter housels posting on comp.os.minix.
* thanks peter!
*
* mjr: 68881 version added
*/
#if !defined (__M68881__) && !defined (sfp004)
#if !(defined(unix) || defined(minix))
#include <stddef.h>
#include <stdlib.h>
#include <float.h>
#endif
#include <errno.h>
#include <assert.h>
#include <math.h>
#ifdef minix
#include "lib.h"
#endif
#ifndef _COMPILER_H
#include <compiler.h>
#endif
extern int errno;
#if (defined(unix) || defined(minix))
#define NULL ((void *)0)
#endif
#define Ise(c) ((c == 'e') || (c == 'E') || (c == 'd') || (c == 'D'))
#define Isdigit(c) ((c <= '9') && (c >= '0'))
#define Isspace(c) ((c == ' ') || (c == '\t'))
#define Issign(c) ((c == '-') || (c == '+'))
#define IsValidTrail(c) ((c == 'F') || (c == 'L'))
#define Val(c) ((c - '0'))
#define MAXDOUBLE DBL_MAX
#define MINDOUBLE DBL_MIN
#define MAXF 1.797693134862316
#define MINF 2.225073858507201
#define MAXE 308
#define MINE (-308)
/* another alias for ieee double */
struct ldouble {
unsigned long hi, lo;
};
static int __ten_mul __PROTO((double *acc, int digit));
static double __adjust __PROTO((double *acc, int dexp, int sign));
#ifdef __OLD__
static double __ten_pow __PROTO((double r, int e));
#endif
/*
* mul 64 bit accumulator by 10 and add digit
* algorithm:
* 10x = 2( 4x + x ) == ( x<<2 + x) << 1
* result = 10x + digit
*/
static int __ten_mul(acc, digit)
double *acc;
int digit;
{
register unsigned long d0, d1, d2, d3;
register short i;
d2 = d0 = ((struct ldouble *)acc)->hi;
d3 = d1 = ((struct ldouble *)acc)->lo;
/* check possibility of overflow */
/* if( (d0 & 0x0FFF0000L) >= 0x0ccc0000L ) */
/* if( (d0 & 0x70000000L) != 0 ) */
if( (d0 & 0xF0000000L) != 0 )
/* report overflow somehow */
return 1;
/* 10acc == 2(4acc + acc) */
for(i = 0; i < 2; i++)
{ /* 4acc == ((acc) << 2) */
asm volatile(" lsll #1,%1;
roxll #1,%0" /* shift L 64 bit acc 1bit */
: "=d" (d0), "=d" (d1)
: "0" (d0), "1" (d1) );
}
/* 4acc + acc */
asm volatile(" addl %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
asm volatile(" addxl %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
/* (4acc + acc) << 1 */
asm volatile(" lsll #1,%1;
roxll #1,%0" /* shift L 64 bit acc 1bit */
: "=d" (d0), "=d" (d1)
: "0" (d0), "1" (d1) );
/* add in digit */
d2 = 0;
d3 = digit;
asm volatile(" addl %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
asm volatile(" addxl %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
/* stuff result back into acc */
((struct ldouble *)acc)->hi = d0;
((struct ldouble *)acc)->lo = d1;
return 0; /* no overflow */
}
#include "flonum.h"
static double __adjust(acc, dexp, sign)
double *acc; /* the 64 bit accumulator */
int dexp; /* decimal exponent */
int sign; /* sign flag */
{
register unsigned long d0, d1, d2, d3;
register short i;
register short bexp = 0; /* binary exponent */
unsigned short tmp[4];
double r;
#ifdef __STDC__
double __normdf( double d, int exp, int sign, int rbits );
double ldexp(double d, int exp);
#else
extern double __normdf();
extern double ldexp();
#endif
d0 = ((struct ldouble *)acc)->hi;
d1 = ((struct ldouble *)acc)->lo;
while(dexp != 0)
{ /* something to do */
if(dexp > 0)
{ /* need to scale up by mul */
while(d0 > 429496729 ) /* 2**31 / 5 */
{ /* possibility of overflow, div/2 */
asm volatile(" lsrl #1,%1;
roxrl #1,%0"
: "=d" (d1), "=d" (d0)
: "0" (d1), "1" (d0));
bexp++;
}
/* acc * 10 == 2(4acc + acc) */
d2 = d0;
d3 = d1;
for(i = 0; i < 2; i++)
{ /* 4acc == ((acc) << 2) */
asm volatile(" lsll #1,%1;
roxll #1,%0" /* shift L 64 bit acc 1bit */
: "=d" (d0), "=d" (d1)
: "0" (d0), "1" (d1) );
}
/* 4acc + acc */
asm volatile(" addl %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
asm volatile(" addxl %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
/* (4acc + acc) << 1 */
bexp++; /* increment exponent to effectively acc * 10 */
dexp--;
}
else /* (dexp < 0) */
{ /* scale down by 10 */
while((d0 & 0xC0000000L) == 0)
{ /* preserve prec by ensuring upper bits are set before div */
asm volatile(" lsll #1,%1;
roxll #1,%0" /* shift L to move bits up */
: "=d" (d0), "=d" (d1)
: "0" (d0), "1" (d1) );
bexp--; /* compensate for the shift */
}
/* acc/10 = acc/5/2 */
*((unsigned long *)&tmp[0]) = d0;
*((unsigned long *)&tmp[2]) = d1;
d2 = (unsigned long)tmp[0];
asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
tmp[0] = (unsigned short)d2; /* the quotient only */
for(i = 1; i < 4; i++)
{
d2 = (d2 & 0xFFFF0000L) | (unsigned long)tmp[i]; /* REM|next */
asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
tmp[i] = (unsigned short)d2;
}
d0 = *((unsigned long *)&tmp[0]);
d1 = *((unsigned long *)&tmp[2]);
/* acc/2 */
bexp--; /* div/2 taken care of by decrementing binary exp */
dexp++;
}
}
/* stuff the result back into acc */
((struct ldouble *)acc)->hi = d0;
((struct ldouble *)acc)->lo = d1;
/* normalize it */
r = __normdf( *acc, ((0x03ff - 1) + 53), (sign)? -1 : 0, 0 );
/* now shove in the binary exponent */
return ldexp(r, bexp);
}
/* flags */
#define SIGN 0x01
#define ESIGN 0x02
#define DECP 0x04
#define CONVF 0x08
double strtod (s, endptr)
const register char *s;
char **endptr;
{
double accum = 0.0;
register short flags = 0;
register short texp = 0;
register short e = 0;
double zero = 0.0;
const char *tmp;
assert ((s != NULL));
if(endptr != NULL) *endptr = (char *)s;
while(Isspace(*s)) s++;
if(*s == '\0')
{ /* just leading spaces */
return zero;
}
if(Issign(*s))
{
if(*s == '-') flags = SIGN;
if(*++s == '\0')
{ /* "+|-" : should be an error ? */
return zero;
}
}
do {
if (Isdigit(*s))
{
flags |= CONVF;
if( __ten_mul(&accum, Val(*s)) ) texp++;
if(flags & DECP) texp--;
}
else if(*s == '.')
{
if (flags & DECP) /* second decimal point */
break;
flags |= DECP;
}
else
break;
s++;
} while (1);
if(Ise(*s))
{
if(*++s != '\0') /* skip e|E|d|D */
{ /* ! ([s]xxx[.[yyy]]e) */
tmp = s;
while(Isspace(*s)) s++; /* Ansi allows spaces after e */
if(*s != '\0')
{ /* ! ([s]xxx[.[yyy]]e[space]) */
if(Issign(*s))
if(*s++ == '-') flags |= ESIGN;
if(*s != '\0')
{ /* ! ([s]xxx[.[yyy]]e[s]) -- error?? */
for(; Isdigit(*s); s++)
e = (((e << 2) + e) << 1) + Val(*s);
if(IsValidTrail(*s)) s++;
/* dont care what comes after this */
if(flags & ESIGN)
texp -= e;
else
texp += e;
}
}
if(s == tmp) s++; /* back up pointer for a trailing e|E|d|D */
}
}
if((endptr != NULL) && (flags & CONVF))
*endptr = (char *) s;
if(accum == zero) return zero;
return __adjust(&accum, (int)texp, (int)(flags & SIGN));
}
double atof(s)
const char *s;
{
#ifdef __OLD__
extern double strtod();
#endif
return strtod(s, (char **)NULL);
}
/* old naive coding */
#ifdef __OLD__
static double __ten_pow(r, e)
double r;
register int e;
{
if(e < 0)
for(; e < 0; e++) r /= 10.0;
else
for(; e > 0; --e) r *= 10.0;
return r;
}
#define RET(X) {goto ret;}
double strtod (s, endptr)
const register char *s;
char **endptr;
{
double f = 0.1, r = 0.0, accum = 0.0;
register int e = 0, esign = 0, sign = 0;
register int texp = 0, countexp;
assert ((s != NULL));
while(Isspace(*s)) s++;
if(*s == '\0') RET(r); /* just spaces */
if(Issign(*s))
{
if(*s == '-') sign = 1;
s++;
if(*s == '\0') RET(r); /* "+|-" : should be an error ? */
}
countexp = 0;
while(Isdigit(*s))
{
if(!countexp && (*s != '0')) countexp = 1;
accum = (accum * 10.0) + Val(*s);
/* should check for overflow here somehow */
s++;
if(countexp) texp++;
}
r = (sign ? (-accum) : accum);
if(*s == '\0') RET(r); /* [s]xxx */
countexp = (texp == 0);
if(*s == '.')
{
s++;
if(*s == '\0') RET(r); /* [s]xxx. */
if(!Ise(*s))
{
while(Isdigit(*s))
{
if(countexp && (*s == '0')) --texp;
if(countexp && (*s != '0')) countexp = 0;
accum = accum + (Val(*s) * f);
f = f / 10.0;
/* overflow (w + f) ? */
s++;
}
r = (sign ? (-accum) : accum);
if(*s == '\0') RET(r); /* [s]xxx.yyy */
}
}
if(!Ise(*s)) RET(r); /* [s]xxx[.[yyy]] */
s++; /* skip e */
if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e */
while(Isspace(*s)) s++;
if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e[space] */
if(Issign(*s))
{
if(*s == '-') esign = 1;
s++;
if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e[s] */
}
while(Isdigit(*s))
{
e = (e * 10) + Val(*s);
s++;
}
/* dont care what comes after this */
if(esign) e = (-e);
texp += e;
/* overflow ? */ /* FIXME */
if( texp > MAXE)
{
if( ! ((texp == (MAXE+1)) && (accum <= MAXF)))
{
errno = ERANGE;
r = ((sign) ? -HUGE_VAL : HUGE_VAL);
RET(r);
}
}
/* underflow -- in reality occurs before this */ /* FIXME */
if(texp < MINE)
{
errno = ERANGE;
r = ((sign) ? -HUGE_VAL : HUGE_VAL);
RET(r);
}
r = __ten_pow(r, e);
/* fall through */
/* all returns end up here, with return value in r */
ret:
if(endptr != NULL)
*endptr = s;
return r;
}
#endif /* __OLD__ */
#else __M68881__ || sfp004
#if 0
#ifndef sfp004
# define _M68881 /* use the predefined inline functions */
#endif sfp004
#endif 0
/* M.R's kludgy atof --- 881 version. */
/* uses long integer accumulators and extended precision to put them */
/* together in the fpu. The conversion long to extended is done completely */
/* on the 881. */
/* using extended precision in _float_ avoids rounding errors. */
/* 12.7.1989, 11.10.90, 28.1.91, 24.11.91 */
/* On overflow, only +-infinity is returned (the 68881's default). */
/* 24.11.91: return +-MAXDOUBLE instead of +- INFINITY or NAN */
/* set errno to ERANGE/EDOM */
# include <ctype.h>
# include <stdio.h>
# include <float.h>
# include <math.h>
# include <errno.h>
# include "flonum.h"
double atof( const char * );
double strtod( const char *, const char ** );
double _Float_( long, long, long, long );
# define true 1
# define false 0
# define CharIsDigit ( isdigit(*Text) )
# define Digit ((*Text-'0'))
#if 0
static unsigned long
__notanumber[2] = { 0x7fffffffL, 0xffffffffL }; /* ieee NAN */
# define NAN (*((double *)&__notanumber[0]))
# endif 0
# define ten_mul(X) ((((X) << 2) + (X)) << 1)
double strtod( const char * Save, const char ** Endptr )
{
register int Count; int Negative = false, ExpNegative = false;
double Value;
union double_di * l_Value;
register long Exponent, Exp_Temp;
register long Value_1, Value_2;
register char c;
register char * Text;
register char * Places;
char Buffer[15];
l_Value = (union double_di *) &Value;
Text = Save;
Places = Buffer;
/* skip over leading whitespace */
while (isspace(*Text)) Text++;
if (*Text == '-') {
Negative = true;
Text++;
} else
if (*Text == '+') {
Negative = false;
Text++;
} else
if( *Text == 0 ) {
if( Endptr != NULL ) *Endptr = Text;
return 0.0;
}
/* Process the 'f'-part */
/* ignore any digit beyond the 15th */
/* BUG: may overflow if more than 10 digits precede the '.' */
/* to cure use to accumulators as is being done for the digits after */
/* the '.'. */
Exp_Temp = 0; /* needed later on for the exponential part */
Value_1 = 0; Value_2 = 0; Count = 0; Exponent = 0;
while( CharIsDigit ) { /* process digits before '.' */
if( Count < 15 ) {
Count++;
*Places++ = Digit;
}
Text++;
}
if ( *Text == '.') {
Text++;
while( CharIsDigit ) { /* process digits after '.' */
if( Count < 15 ) {
Count++;
*Places++ = Digit;
Exponent--;
}
Text++;
}
}
Places = Buffer;
/* Now, Places points to a vector of <= 15 digits */
/* text points to the position immediately after the end of the mantissa */
/* Value_2 will contain the equiv. of the 8 least significant digits, while */
/* Value_1 will contain the equiv. of the 7 most significant digits (if any) */
/* and therefore has to be multiplied by 10^8 */
/* no overflow possible in the temporary buffers */
while( Count > 8 ) {
Value_1 = ten_mul( Value_1 ); Value_1 += *Places++;
Count--;
}
while( Count > 0 ) {
Value_2 = ten_mul( Value_2 ); Value_2 += *Places++;
Count--;
}
/* 'e'-Part */
if ( *Text == 'e' || *Text == 'E' || *Text == 'd' || *Text == 'D' ) {
char * Tail = Text;
Text++;
/* skip over whitespace since ANSI allows space after e|E|d|D */
while (isspace(*Text)) Text++;
if ( * Text == '-' ) {
ExpNegative = true;
Text++;
} else
if( * Text == '+' ) {
ExpNegative = false;
Text++;
}
if( !CharIsDigit ) {
*Endptr = Tail; /* take the 'e|E|d|D' as part of the characters */
goto Ende; /* following the number */
} else {
/* Exponents may have at most 3 digits, everything beyond this will be */
/* silently ignored */
Count = 0;
while( CharIsDigit && (Count < 3) ) {
Exp_Temp = ten_mul( Exp_Temp ); Exp_Temp += Digit;
Count++;
Text++;
}
if( ExpNegative ) Exp_Temp = -Exp_Temp;
Exponent += Exp_Temp;
}
}
Value = _Float_( Value_1, Exponent+8L, Value_2, Exponent );
if( Endptr != NULL ) *Endptr = Text;
if( (l_Value -> i[0]) >= 0x7ff00000U ) { /* INFINITY or NAN */
fprintf(stderr," strtod: OVERFLOW error\n");
errno = ERANGE;
Value = DBL_MAX;
}
Ende:
if( Negative ) {
Value = -Value;
}
return( Value );
# if 0
Error:
fputs("\njunk number \"",stderr); fputs(Save,stderr);
fputs("\" --- returning NAN\n",stderr);
errno = ERANGE;
if( Endptr != NULL ) *Endptr = Text;
return(NAN); /* == Not A Number (NAN) */
# endif
}
double atof( const char * Text )
{
return(strtod(Text,(char **)NULL));
}
/*
* double _Float_( long Value_1, long Exponent_1, long Value_2, long Exponent_2 )
*
* merges the accumulators Value_1, Value_2 and the Exponent to a double
* precision float
* called by strtod()
*
* does all floating point computations with extended precision on the fpu
*/
#endif __M68881__ || sfp004
#ifdef __M68881__
asm(
".even
.text
.globl __Float_
__Float_:
ftentoxl a7@(8),fp1 | load Exponent_1
ftentoxl a7@(16),fp0 | load Exponent_2
fmull a7@(12),fp0 | fmull Value_2 -> fp0
fmull a7@(4),fp1 | fmull Value_1 -> fp1
faddx fp0,fp1
fmoved fp1,a7@- | fmoved fp1 -> d0/d1
moveml a7@+,d0-d1
rts
");
#endif __M68881
#ifdef sfp004
asm("| mjr, 30.1.1991
|
| base = 0xfffa50
| the fpu addresses are taken relativ to 'base':
|
| a0: fpu base address
|
| waiting loop ...
|
| wait:
| ww: cmpiw #0x8900,a0@(resp)
| beq ww
| is coded directly by
| .long 0x0c688900, 0xfff067f8 (a0)
| and
| www: tst.w a0@(resp)
| bmi.b www
| is coded by
| .word 0x4a68,0xfff0,0x6bfa | test
|
comm = -6
resp = -16
zahl = 0
.globl __Float_
.even
.text
__Float_:
lea 0xfffffa50:w,a0 | fpu address
movew #0x4092,a0@(comm) | ftentoxl -> fp1
.long 0x0c688900, 0xfff067f8
movel a7@(8),a0@ | load Exponent_1
movew #0x4112,a0@(comm) | ftentoxl -> fp2
.long 0x0c688900, 0xfff067f8
movel a7@(16),a0@ | load Exponent_2
movew #0x4123,a0@(comm) | fmull Value_2 -> fp2
.long 0x0c688900, 0xfff067f8
movel a7@(12),a0@ | load Value_2
movew #0x40a3,a0@(comm) | fmull Value_1 -> fp1
.long 0x0c688900, 0xfff067f8
movel a7@(4),a0@ | load Value_1
movew #0x08a2,a0@(comm) | faddx fp2 -> fp1
.word 0x4a68,0xfff0,0x6bfa | test
movew #0x7480,a0@(comm) | fmoved fp1 -> d0/d1
.long 0x0c688900, 0xfff067f8
movel a0@,d0
movel a0@,d1
rts
");
#endif sfp004
#ifdef TEST
#if 0
#ifdef __MSHORT__
#error "please run this test in 32 bit int mode"
#endif
#endif
#define NTEST 10000L
#ifdef __MSHORT__
# define RAND_MAX (0x7FFF) /* maximum value from rand() */
#else
# define RAND_MAX (0x7FFFFFFFL) /* maximum value from rand() */
#endif
main()
{
double expected, result, e, max_abs_err;
char buf[128];
register long i, errs;
register long s;
#ifdef __STDC__
double atof(const char *);
int rand(void);
#else
extern double atof();
extern int rand();
#endif
#if 0
expected = atof("3.14159265358979e23");
expected = atof("3.141");
expected = atof(".31415");
printf("%f\n\n", expected);
expected = atof("3.1415");
printf("%f\n\n", expected);
expected = atof("31.415");
printf("%f\n\n", expected);
expected = atof("314.15");
printf("%f\n\n", expected);
expected = atof(".31415");
printf("%f\n\n", expected);
expected = atof(".031415");
printf("%f\n\n", expected);
expected = atof(".0031415");
printf("%f\n\n", expected);
expected = atof(".00031415");
printf("%f\n\n", expected);
expected = atof(".000031415");
printf("%f\n\n", expected);
expected = atof("-3.1415e-9");
printf("%20.15e\n\n", expected);
expected = atof("+3.1415e+009");
printf("%20.15e\n\n", expected);
#endif
expected = atof("+3.123456789123456789");
printf("%30.25e\n\n", expected);
expected = atof(".000003123456789123456789");
printf("%30.25e\n\n", expected);
expected = atof("3.1234567891234567890000000000");
printf("%30.25e\n\n", expected);
expected = atof("9.22337999999999999999999999999999999999999999");
printf("%47.45e\n\n", expected);
expected = atof("1.0000000000000000000");
printf("%25.19e\n\n", expected);
expected = atof("1.00000000000000000000");
printf("%26.20e\n\n", expected);
expected = atof("1.000000000000000000000");
printf("%27.21e\n\n", expected);
expected = atof("1.000000000000000000000000");
printf("%30.24e\n\n", expected);
#if 0
expected = atof("1.7e+308");
if(errno != 0)
{
printf("%d\n", errno);
}
else printf("1.7e308 OK %g\n", expected);
expected = atof("1.797693e308"); /* anything gt looses */
if(errno != 0)
{
printf("%d\n", errno);
}
else printf("Max OK %g\n", expected);
expected = atof("2.225073858507201E-307");
if(errno != 0)
{
printf("%d\n", errno, expected);
}
else printf("Min OK %g\n", expected);
#endif
max_abs_err = 0.0;
for(errs = 0, i = 0; i < NTEST; i++)
{
expected = (double)(s = rand()) / (double)rand();
if(s > (RAND_MAX >> 1)) expected = -expected;
sprintf(buf, "%.14e", expected);
result = atof(buf);
e = (expected == 0.0) ? result : (result - expected)/expected;
if(e < 0) e = (-e);
if(e > 1.0e-6)
{
errs++; printf("%.14e %s %.14e (%.14e)\n", expected, buf, result, e);
}
if (e > max_abs_err) max_abs_err = e;
}
printf("%ld Error(s), Max abs err %.14e\n", errs, max_abs_err);
}
#endif /* TEST */